Add lambert azimuthal equal area coord_system#1738
Conversation
lib/iris/coord_systems.py
Outdated
There was a problem hiding this comment.
Please could you add "Defaults to 0." for these arguments.
|
👍 Now it just needs a new Cartopy release... |
|
It probably also needs an entry to the what's new. |
Yes, good point. 👍 |
951b4a3 to
5441678
Compare
|
@bjlittle , this is the pull request I was talking about that is now needed. |
lib/iris/fileformats/netcdf.py
Outdated
| (iris.coord_systems.TransverseMercator, | ||
| iris.coord_systems.LambertAzimuthalEqualArea)): | ||
| units = 'm' | ||
|
|
There was a problem hiding this comment.
This block was removed in PR #1943 but I don't agree that it should have been removed.
Section 4.1 and 4.2 of the cf conventions says:
Coordinates of latitude with respect to a rotated pole should be given units of degrees
also 'm' is a perfectly fine unit (see the example in this section).
There was a problem hiding this comment.
@lbdreyer Interesting that swapping this back doesn't cause any test failures ...
There was a problem hiding this comment.
In fact, if you look at Example 5.6, that's a rotated_lat_lon grid with units of "degrees". Also Example 5.10 is a tranverse_mercator grid with units of "m".
There was a problem hiding this comment.
Right, I think I must have gotten confused trying to rebase a PR that was over a year old!
On closer inspection I now see that, in the case of tranverse_mercator we're just passing the units from the coord. Although this has it's disadvantages (i.e. it means that a user could produce a netcdf file with a badly defined coordinate:
e.g. DimCoord(array([30, 60]), standard_name='latitude', units=Unit('m'), coord_system=RotatedGeogCS(10.0, 20.0))
But I would prefer for users to have this flexibility, and cf says that there is no default value for lat or lon units.
As the CF conventions recommend a specific unit for GeogCS and RotatedGeogCS I think you could argue that we should 'enforce' these:
The recommended unit of latitude is degrees_north
and
Coordinates of latitude with respect to a rotated pole should be given units of degrees
|
Failing because of this: |
|
That failure was fixed in #2144, targeted at the v1.10.x branch. This change needs to make its way onto master, and then you can rebase onto that to get passing tests here (or we can trust that this feature doesn't cause the failure). There is another test failure too (an error actually), not sure what that one is caused by. |
|
I am going to add a test for the coordinate units discussed in the above comment (I'm hoping that's not a controversial change), and then rebase. |
5441678 to
0af35ac
Compare
| def check(self, coord_system, units, expected_units): | ||
| coord = iris.coords.DimCoord([30, 45], 'latitude', units=units, | ||
| coord_system=coord_system) | ||
| saver = Saver('filename', 'NETCDF4') |
There was a problem hiding this comment.
Argh! This is wrong. I was going to mock it out.
I had made _cf_coord_identity a static method in the hopes that I wouldn't have to instantiate the Saver class to test it, but I was getting an error
>>> coord = iris.coords.DimCoord([45], 'latitude', units='m')
>>> r = Saver._cf_coord_identity(coord)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: unbound method _cf_coord_identity() must be called with Saver instance as first argument (got DimCoord instance instead)
35883f7 to
9e3a3b8
Compare
|
@bjlittle The tests are now passing and the static method is working correctly. |
|
Hi @lbdreyer i'm afraid something has occurred on master such that your branch now has conflicts please rebase and fix so we can see a plausible change with test results are you still waiting on a cartopy release for this to be able to work properly? mark |
9e3a3b8 to
cde0ac4
Compare
|
I have rebased and the tests are passing again.
No, the change that this required was in the Cartopy 0.14 release |
bjlittle
left a comment
There was a problem hiding this comment.
Some minor changes, and some feedback from @djkirkham and @marqh would be welcome regarding rotated pole and (CF spec) units of degrees
lib/iris/fileformats/netcdf.py
Outdated
| (iris.coord_systems.TransverseMercator, | ||
| iris.coord_systems.LambertAzimuthalEqualArea)): | ||
| units = 'm' | ||
|
|
There was a problem hiding this comment.
@lbdreyer Interesting that swapping this back doesn't cause any test failures ...
lib/iris/fileformats/netcdf.py
Outdated
|
|
||
| return (coord.standard_name, coord.long_name, units) | ||
| elif isinstance(coord.coord_system, iris.coord_systems.RotatedGeogCS): | ||
| units = 'degrees' |
There was a problem hiding this comment.
@lbdreyer As per the CF spec sections 4.1 and 4.2.
@djkirkham and @marqh This reverts your changes from #1943. Can we all agree this is what we want for rotated pole units, as per CF.
There was a problem hiding this comment.
If I remember correctly, we took this out because the units may be radians, say, so changing to degrees would be incorrect.
There was a problem hiding this comment.
Thanks @djkirkham, that makes sense. I will remove this.
| self.semi_major_axis = 6377563.396 | ||
| self.semi_minor_axis = 6356256.909 | ||
| self.ellipsoid = GeogCS(self.semi_major_axis, self.semi_minor_axis) | ||
| self.laea_cs = LambertAzimuthalEqualArea( |
There was a problem hiding this comment.
@lbdreyer Add in the false_easting and false_northing ... thanks!
| def check_call(self, coord_system, units, expected_units): | ||
| coord = iris.coords.DimCoord([30, 45], 'latitude', units=units, | ||
| coord_system=coord_system) | ||
| result = Saver._cf_coord_identity(coord) |
There was a problem hiding this comment.
@lbdreyer Hence the need for @staticmethod, okay ... 😄
| class Test__cf_coord_identity(tests.IrisTest): | ||
| def check_call(self, coord_system, units, expected_units): | ||
| coord = iris.coords.DimCoord([30, 45], 'latitude', units=units, | ||
| coord_system=coord_system) |
There was a problem hiding this comment.
@lbdreyer Ideally, to give full coverage, could you also test longitude also ... thanks!
| self.check_call(coord_system=crs, units='degrees', | ||
| expected_units='degrees') | ||
|
|
||
| def test_crs_with_no_default_units(self): |
There was a problem hiding this comment.
@lbdreyer This is really a pass-thru rather than default behaviour ... what do you think?
| latitude_of_projection_origin=52, | ||
| longitude_of_projection_origin=10, | ||
| false_easting=100, | ||
| false_northing=200) |
There was a problem hiding this comment.
@lbdreyer Did you want to pass-thru a choice of ellipsoid here ...
5bbe26e to
fa72ee7
Compare
|
@lbdreyer Great work! Thanks 👍 |
Depends on SciTools/iris-test-data#37 and SciTools/cartopy#619